### How Innovation in Political Participation Could Increase Perceived Legitimacy
### Replication Code

### Leopold Weil, Regula Hänggli
### 17 May 2021

###################

## Initialize libraries

library(readr)
library(magrittr)
library(tidyverse)
library(stargazer)
library(estimatr)
library(haven)
library(cjoint)
library(sandwich)
library(psych)

## Load data

data1 <- read_csv("replication_data.csv")


## Construct indices: political trust, general trust, political efficacy

data1 %<>%
  mutate(efficacy =  efficacy_int_1 + efficacy_int_2 - efficacy_int_3 + efficacy_ext_1 + efficacy_ext_2 - efficacy_ext_3) %>%
  mutate(pol.trust = pol.trust1 + pol.trust2 + pol.trust3 + pol.trust4 + pol.trust5 + pol.trust6 + pol.trust7 + pol.trust8 + pol.trust9) %>%
  mutate(gen.trust = general.trust1 + general.trust2 + general.trust3)


## Randomization checks: ANOVA

gender.table <- table(data1$gender, data1$treatment)
fisher.test(gender.table, simulate.p.value = T)

random.check.vars <- list(data1$age, data1$pol.interest, data1$leftright, data1$pol.trust, data1$gen.trust)

for (i in 1:length(random.check.vars)){
  aov(random.check.vars[[i]]~data1$treatment) %>%
    summary %>%
    print
}


## Descriptive sample plots

ggplot(data = data1, aes(x = gender)) +
  geom_bar(aes(y = stat(count/3))) +
  coord_flip()+
  labs(y = "Number of participants", x = NULL) +
  theme_bw()

ggplot(data = data1, aes(x = age)) +
  geom_histogram(aes(y = stat(count/3)), color = "white")+
  labs(x = "Age", y = "Number of participants") +
  theme_bw()

ggplot(data = data1, aes(x = educ)) +
  geom_bar(aes(y = stat(count/3))) +
  coord_flip() +
  labs(y = "Number of participants", x = NULL) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust = 0))


## PCA of legitimacy items

data2 <- data1 %>%
  dplyr::select(acceptance01, acceptance02, fairness01, fairness02, influence01, influence02, trust01, trust02)

data2 <- data2 - 1

stargazer(as.data.frame(data2), type = "text", header = F, summary.stat = c("mean", "sd", "median", "max","min","n"))

pca1 <- pca(data2, nfactors = 8, rotate = "none")

pc1 <- data.frame(PC1 = pca1$loadings[,1])
stargazer(pc1, summary = F, header = F, type="text")

pc_all <- data.frame(PC1 = pca1$loadings[,1],
                     PC2 = pca1$loadings[,2],
                     PC3 = pca1$loadings[,3],
                     PC4 = pca1$loadings[,4],
                     PC5 = pca1$loadings[,5],
                     PC6 = pca1$loadings[,6],
                     PC7 = pca1$loadings[,7],
                     PC8 = pca1$loadings[,8]
)

stargazer(pc_all, summary = F, header = F, type="text")

pc_all2 <- pca1$Vaccounted[1:3,] %>%
  as.data.frame
row.names(pc_all2)[1] <- "Eigenvalue"

stargazer(pc_all2, summary = F, header = F, type="text")

data1$legitimacy <- pca1$scores[,1]


## Conjoint models

# Baseline model

data3 <- data1 %>%
  dplyr::select(legitimacy, agendasetting.fac, voting, outcome.fac, ID) %>%
  na.omit

cj1 <- amce(legitimacy ~ agendasetting.fac + voting +  outcome.fac, data = data3, respondent.id = "ID")

summary(cj1)

# Baseline model with different baseline category for voting

voting.baseline = list(voting = "majority")

cj1.1 <- amce(legitimacy ~ agendasetting.fac + voting +  outcome.fac, data = data3, respondent.id = "ID", baselines = voting.baseline)

summary(cj1.1)

## Plots for conjoint models

# Initialize data frame for plotting
cj.data <- data.frame(label = c("Voting:", "   majority", "   preferential", "   quadratic", "Agendasetting:", "   bottom-up", "Outcome:", "   preferred"), AMCE = rep(NA,8), se = rep(NA,8))

cj.data$label <- factor(cj.data$label, levels = rev(cj.data$label))

cj.data[6,2:3] <- cj1$estimates$agendasettingfac
cj.data[2:4,2:3] <- t(cj1$estimates$voting)
cj.data[8, 2:3] <- cj1$estimates$outcomefac

# Define plotting function
conjoint.plot <- function(data){
  ggplot(data = data) +
    geom_pointrange(aes(x=label, y=AMCE, ymin = AMCE-1.96*se, ymax = AMCE+1.96*se), color = "steelblue") +
    coord_flip() +
    theme_bw() +
    theme(axis.text.y = element_text(hjust = 0)) +
    labs(x = NULL, y = "Average Marginal Component Effect") +
    geom_hline(yintercept = 0, color = "firebrick", linetype = "dashed")
}

# Generate plot
conjoint.plot(cj.data)

# Initialize second data frame
cj.data2 <- cj.data
cj.data2[2:4,2:3] <- t(cj1.1$estimates$voting)
cj.data2$label <- factor(c("Voting:", "   none", "   preferential", "   quadratic", "Agendasetting:", "   bottom-up", "Outcome:", "   preferred"), levels = rev(c("Voting:", "   none", "   preferential", "   quadratic", "Agendasetting:", "   bottom-up", "Outcome:", "   preferred")))

# Generate plot
conjoint.plot(cj.data2)


## Conjoint models for subgroups (interaction)

# Estimate models
cj.int.1 <- amce(legitimacy ~ agendasetting.fac + voting, data = data3, respondent.id = "ID", subset = data3$outcome.fac == "preferred")

summary(cj.int.1)

cj.int.2 <- amce(legitimacy ~ agendasetting.fac + voting, data = data3, respondent.id = "ID", subset = data3$outcome.fac == "opposed")

summary(cj.int.2)

# Initialize plot data frame
cj.data3 <- data.frame(label = rep(c("Voting:", "   majoritarian", "   preferential", "   quadratic", "Agendasetting:", "   bottom-up"), 2), AMCE = rep(NA,12), se = rep(NA,12), outcome = c(rep("opposed", 6), rep("preferred", 6)))

cj.data3$label <- factor(cj.data3$label, levels = rev(c("Voting:", "   majoritarian", "   preferential", "   quadratic", "Agendasetting:", "   bottom-up")))

cj.data3[6,2:3] <- cj.int.2$estimates$agendasettingfac
cj.data3[12,2:3] <- cj.int.1$estimates$agendasettingfac

cj.data3[2:4,2:3] <- t(cj.int.2$estimates$voting)
cj.data3[8:10,2:3] <- t(cj.int.1$estimates$voting)

# Generate plot
ggplot(data = cj.data3) +
  geom_pointrange(aes(x=label, y=AMCE, ymin = AMCE-1.96*se, ymax = AMCE+1.96*se, color = outcome), position = position_dodge(width = 0.3)) +
  coord_flip() +
  theme_bw() +
  theme(axis.text.y = element_text(hjust = 0)) +
  labs(x = NULL, y = "Average Marginal Component Effect conditional on outcome favourability") +
  geom_hline(yintercept = 0, color = "firebrick", linetype = "dashed") 

## Predicted value computations

# Estimate model and clustered standard errors
lin.mod <- lm(legitimacy ~ (agendasetting.fac + voting) * outcome.fac, data = data1)
vcov.clst <- vcovHC(lin.mod,type="HC3")

# Define computation function
simulation <- function(covariates){
  
  S <- mvrnorm(1000, coef(lin.mod), vcov.clst)
  
  x <- matrix(covariates)
  
  mu <- S %*% x
  
  quantile(mu, c(0.025, 0.5, 0.975)) %>%
    return
  
}

# Define covariate scenarios
covariates0 <- c(1, #intercept
                 rep(0, 9)) 


covariates1 <- c(1, #intercept
                 1, #agendasetting bottom up
                 0, # majority vote
                 0, # preferential vote
                 0, # quadratic vote
                 0, # outcome preferred
                 0, # agendasetting bottum up x outcome preferred
                 0, # majority vote x outcome preferred
                 0, # preferential vote x outcome preferred
                 0) # quadratic vote x outcome preferred

covariates2 <- c(1, #intercept
                 0, #agendasetting bottom up
                 0, # majority vote
                 1, # preferential vote
                 0, # quadratic vote
                 0, # outcome preferred
                 0, # agendasetting bottum up x outcome preferred
                 0, # majority vote x outcome preferred
                 0, # preferential vote x outcome preferred
                 0) # quadratic vote x outcome preferred

covariates3 <- c(1, #intercept
                 1, #agendasetting bottom up
                 0, # majority vote
                 1, # preferential vote
                 0, # quadratic vote
                 1, # outcome preferred
                 0, # agendasetting bottum up x outcome preferred
                 0, # majority vote x outcome preferred
                 1, # preferential vote x outcome preferred
                 0) # quadratic vote x outcome preferred


# Run computations and attach results
sims <- simulation(covariates0)
sims <- rbind(sims, simulation(covariates1))
sims <- rbind(sims, simulation(covariates2))
sims <- rbind(sims, simulation(covariates3)) %>%
  as.data.frame

# Plot results
scenarios <- c("Baseline", "1", "2", "3")
sims$szenario <- factor(scenarios, levels = scenarios)

ggplot(data = sims) +
  geom_pointrange(aes(x = szenario, y = `50%`, ymin = `2.5%`, ymax= `97.5%`), color = "steelblue") +
  labs(y = "Legitimacy index", x = "Scenario") +
  theme_bw()